function varargout = slscatter(X, g, type, varargin)
%SLSCATTER Compute the scatter matrix 
%
% [ Syntax ]
%   - Sw = slscatter(X, nums, 'Sw', ...)
%   - Sb = slscatter(X, nums, 'Sb', ...)
%   - St = slscatter(X, nums, 'St', ...)
%   - [Sb, Sw] = slscatter(X, nums, 'Sb-Sw', ...)
%   - [St, Sw] = slscatter(X, nums, 'St-Sw', ...)
%   - [St, Sb, Sw] = slscatter(X, nums, 'St-Sb-Sw', ...)
%   - ... = slscatter(X, ginds, ...)
%
% [ Arguments ]
%   - X:        the sample matrix with each column representing a sample
%   - nums:     the number of samples in each classes
%
%   - Sw:       the within-class scattering matrix
%   - Sb:       the between-class scattering matrix
%   - St:       the total scattering matrix
%
% [ Description ]
%   - Sw = slscatter(X, nums, 'Sw', ...) computes the within-class
%     scattering matrix.
%
%     Suppose there are n samples in a d-dimensional space from nc classes,
%     X should be a d x n sample matrix. Please note that the samples
%     from the same class should be grouped together here. nums should be
%     a vector of numbers of length nc specifying the numbers of samples in
%     all the classes.     
%
%   - Sb = slscatter(X, nums, 'Sb', ...) computes the between-class
%     scattering matrix.
%
%   - St = slscatter(X, nums, 'St', ...) computes the total scattering
%     matrix.
%
%   - [Sb, Sw] = slscatter(X, nums, 'Sb-Sw', ...) computes both the
%     between-class scattering matrix and the within-class scattering
%     matrix.
%
%   - [St, Sw] = slscatter(X, nums, 'St-Sw', ...) computes both the total
%     scattering matrix and the within-class scattering matrix.
%
%   - [St, Sb, Sw] = slscatter(X, nums, 'St-Sb-Sw', ...) computes the total
%     scattering matrix, the between-class scattering matrix, and the
%     within-class scattering matrix.
%
%   - ... = slscatter(X, ginds, ...) uses ginds, a cell array of index
%     arrays, to specify the class-ownership of the samples. Here, ginds{k}
%     is an array of indices of the samples from the k-th class. 
%
%     With this syntax, there is no requirement that the samples from the
%     same class are arranged together.
%
%     For all the syntax above, your may specify the following options
%     \{:
%        - sweights:        the vector of sample weights.
%                           By default, it is [], which implies that all
%                           samples share the same weight 1.
%
%        - cmeans:          the class-specific mean vectors, in form of
%                           a d x nc matrix.
%                           By default, it is [], then the means will be
%                           computed from the samples.
%
%        - tmean:           the total mean vector, in form of a d x 1
%                           column vector.
%
%                           By default, it is [], then the total mean will 
%                           be computed from class-specific means.        
%     \:}
%
% [ Remarks ]
%   - The syntax as
%       $ [Sb, Sw] = slscatter(X, nums, 'Sb-Sw', ...); $
%     is functionally equivalent to
%     \{
%          Sb = slscatter(X, nums, 'Sb');
%          Sw = slscatter(X, nums, 'Sw'); 
%     \}
%     However, the former syntax will be more efficient. With the 'Sb-Sw'
%     option, some computation shared by 'Sb' and 'Sw' will be done once,
%     while in the latter syntax they will be done twice.
%     The same rules are also applicable to 'St-Sw' and 'St-Sb-Sw'.
%
% [ History ]
%   - Created by Dahua Lin on Apr 27, 2005
%   - Modified by Dahua Lin on Jul 22, 2007
%       - Enhance the support to standard method
%           - remove the 'pw' method, and thus the method option
%           - remove the dwrule option
%           - add the centers option for pre-computed class centers
%       - Refactor the program structure
%       - Support more flexible syntax
%       - Improve the computation efficiency 
%       
%

%% parse and verify input arguments

error(nargchk(3, inf, nargin));

% X
assert(isnumeric(X) && ndims(X) == 2, 'sltoolbox:slscatter:invalidarg', ...
    'X should be a 2D numeric matrix.');
n = size(X, 2);

% g (nums | ginds)
assert((isnumeric(g) || iscell(g)) && isvector(g), 'sltoolbox:slscatter:invalidarg', ...
    'The 2nd argument can be a numeric vector as nums or a cell vector as ginds.');
nc = length(g);
if ~isequal(size(g), [nc 1])
    g = reshape(g, [nc 1]);
end
if isnumeric(g)
    assert(sum(g) == n, 'sltoolbox:slscatter:nummismatch', ...
        'The number of samples given in nums is not equal to that given in X.');    
    eps = cumsum(g);
    sps = [1; eps(1:end-1)+1];
    ginds = arrayfun(@(sp, ep) sp:ep, sps, eps, 'UniformOutput', false);
else
    ginds = g;
end

% type
assert(ischar(type), 'sltoolbox:slscatter:invalidarg', ...
    'The 3rd argument should be a string indicating what type(s) of scatter matrices are computed.');

% options
opts = struct( ...
    'sweights', [], ...
    'cmeans', [], ...
    'tmean', []);
opts = setopts(opts, varargin{:});

sweights = opts.sweights;
if ~isempty(sweights)
    assert(isnumeric(sweights) && isvector(sweights) && length(sweights) == n, ...
        'sltoolbox:slscatter:invalidopt', 'The option sweights should be a numeric vector of length n.');
    sweights = sweights(:);
end

cmeans = opts.cmeans;
if ~isempty(cmeans)
    assert(isnumeric(cmeans) && isequal(size(cmeans), [d nc]), ...
        'sltoolbox:slscatter:invalidopt', 'The option cmeans should be a 2D matrix of size d x nc.');
end

tmean = opts.tmean;
if ~isempty(tmean)
    assert(isnumeric(tmean) && isequal(size(tmean), [d 1]), ...
        'sltoolbox:slscatter:invalidopt', 'The option tmean should be a vector of size d x 1.');
end


%% main

switch type
    case 'Sw'
        if isempty(cmeans)
            cmeans = compute_cmeans(X, ginds, sweights);
        end
        Sw = compute_Sw(X, ginds, cmeans, sweights);
        varargout = {Sw};
        
    case 'Sb'
        if isempty(cmeans)
            cmeans = compute_cmeans(X, ginds, sweights);
        end
        cweights = compute_cweights(ginds, sweights);
        Sb = slcov(cmeans, cweights, tmean);
        varargout = {Sb};
        
    case 'St'
        if isempty(cmeans)
            St = slcov(X, sweights, tmean);
        else
            Sw = compute_Sw(X, ginds, cmeans, sweights);
            cweights = compute_cweights(ginds, sweights);
            Sb = slcov(cmeans, cweights, tmean);
            St = Sw + Sb;
        end
        varargout = {St};
        
    case 'Sb-Sw'
        if isempty(cmeans)
            cmeans = compute_cmeans(X, ginds, sweights);
        end
        cweights = compute_cweights(ginds, sweights);
        Sw = compute_Sw(X, ginds, cmeans, sweights);
        Sb = slcov(cmeans, cweights, tmean);
        varargout = {Sb, Sw};
        
    case 'St-Sw'
        if isempty(cmeans)
            cmeans = compute_cmeans(X, ginds, sweights);
        end
        cweights = compute_cweights(ginds, sweights);
        Sw = compute_Sw(X, ginds, cmeans, sweights);
        Sb = slcov(cmeans, cweights, tmean);
        St = Sb + Sw;
        varargout = {St, Sw};
        
    case 'St-Sb-Sw'
        if isempty(cmeans)
            cmeans = compute_cmeans(X, ginds, sweights);
        end
        cweights = compute_cweights(ginds, sweights);
        Sw = compute_Sw(X, ginds, cmeans, sweights);
        Sb = slcov(cmeans, cweights, tmean);
        St = Sb + Sw;
        varargout = {St, Sb, Sw};
        
    otherwise
        error('sltoolbox:slscatter:unknowntype', ...
            'The type of scatter matrix named %s is unknown', type);
end


%% computational routines    

function cmeans = compute_cmeans(X, ginds, sweights)

if isempty(sweights)
    cmeans = cellfun(@(u) slmean(X(:, u)), ginds, 'UniformOutput', false);
else
    cmeans = cellfun(@(u) slmean(X(:, u), sweights(u)), ginds, 'UniformOutput', false);
end
cmeans = [cmeans{:}];


function cweights = compute_cweights(ginds, sweights)

if isempty(sweights)
    cweights = cellfun(@numel, ginds);
else
    cweights = cellfun(@(u) sum(sweights(u(:))), ginds);
end


function Sw = compute_Sw(X, ginds, cmeans, sweights)

[D, w] = make_withinclass_diffvecs(X, ginds, cmeans, sweights);
Sw = slcov(D, w, 0);



